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CHAPTER  1 
Introduction  -  Qualitative  Review 

1.1  Motivation 

Interest  in  the  strength  of  naturally  occuring  ice  and 
the  ultimate  size  to  which  marine  structures  must  he  designed 
has  become  a  major  topic  as  the  last  boundaries  of  easily 
accessible  petroleum  resources  are  used.   The  extension  of 
engineering  practice  into  the  world's  colder  regions  is 
subsequently  an  effort  that  is  intimately  related  to  this 
perceived  need. 

To  satisfy  a  curiosity  regarding  this  area  of  offshore 
design,  the  author  set  forth  to  discover  a  rational  approach 
to  establishing  a  meaningful  standard.   It  was  hoped  that 
an  appropriate  theory,  perhaps  based  on  statistical  patterns, 
was  available  that  was  similar  to  the  many  comprehensive 
theories  regarding  wind,  wave,  and  seismic  loading.   Unfor- 
tunately, no  comprehensive  or  universally  accepted  theory 
was  found.  However,  many  prominant  investigators,  including 
Assur  [1  ][2][3],  Bercha  [4]  [5]  [6],  Hirayama  et  al.  [14] 
[15],  and  Korzhavin  [19]  [20]*  have  proposed  numerous  models 
and  approaches  to  the  calculation  of  ice  strength  and  the 
quantitative  prediction  of  necessary  design  size.   Some  of 
these  models  will  be  reviewed  in  subsequent  chapters. 

Consequently,  the  motivation  for  this  thesis  is  to 
recognize  and  segregate  an  accurate  and  useful  theory  for 
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the  prediction  of  ice  pressures. 

1.2  Scope 

The  purpose  of  this  investigation  is  to  analyze  the 
ice-structural  interaction  models  presently  proposed  for  the 
case  of  a  vertical  circular  pile.   The  circular  pile  was 
chosen  primarily  because  it  is  the  most  useful  geometric 
form  in  areas  of  high  tidal  variance.   Other  authors, 
notably  Tryde  [30]  and  Bercha  [5]*  have  examined  the  other 
promising  forms,  the  inclined  wedge  and  the  cone. 

The  investigation  will  take  the  form  of  the  following 
steps : 

1)  Qualitatively  compare  the  available  analytical 
models. 

2)  Either  improve  upon  the  most  promising  model  or 
construct  a  new  model. 

3)  Analyze  sensitivity  to  some  important  parameters. 

4)  Analyze  some  important  ice  properties  affecting 
the  circular  pile  interaction  problem. 

1.3  Structural  Failure 

The  earliest  designs  contemplating  an  ice  environment 
were  highway  and  railroad  bridges  in  North  America  and  the 
Soviet  Union.   Some  authors,  notably  Korzhavin  [19]  and 
Watts  [32],  have  summarized  structural  collapses  known  to 
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be  due  to  ice  forces.   These  failures  have  a  common  thread; 
they  were  all  caused  by  unusual  ice  jams  which,  coupled 
with  an  early  thaw,  produced  unplanned  loadings.   Insights 
gained  from  a  careful  reading  of  these  descriptions  indicate: 

-  very  few  failures  actually  occurred.   This  evidence 
implies  that  early  engineers  had  either  very  accurate 
analytical  skills  (hardly  likely)  or  that  early 
designs  were  overly  conservative. 

-  some  failures  were  by  gradual  ice-erosion. 

-  when  failure  did  occur,  it  was  due  to  an  unusual 
meteorological  complication  or  no  allowance  for 
ice  at  all. 

No  major  marine  structural  failures,  apart  from 
ice-bound  ships,  were  found  in  the  literature.  Perhaps  this  is 
because  this  hostile  environment  has  yet  to  be  seriously 
challenged.   Ice  damage  on  coastline  facilities,  breakwaters, 
docking  facilities,  and  lights,  have  been  recorded.  [26] 
1.4  Application 

Applications  of  ice  engineering  extend  not  only  to  the 
offshore  industry  but  over  a  wide  range  of  related  problems. 
The  knowledge  gained  through  the  solution  of  marine  ice 
uncertainty  can  apply  to  the  following  range  of  actual 
proposals : 

-  over  ice  transportation 

-  ice  airfields 

-  submarine  surfacing  through  ice 
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-  icebreaker  technology  -  navigation 

-  ice  island  construction  for  use  as  a  moblie  platform 

-  iceberg  towing 

The  offshore  oil  industry  has  perhaps  the  greatest 
motivation  in  subduing  arctic  conditions.  Although 
experience  to  date  is  not  extensive,  the  few  structures 
presently  deployed  emphasize  the  continuing  need  for  more 
accurate  analysis. 

Present  North  American  interest  is  concentrated  in 
Cook  Inlet,  Alaska,  [2Q]t   the  North  Slope,  and  the  Beaufort 
Sea  [9]«  Arctic  conditions  are  most  severe,  as  it  will 
be  shown  in  the  next  few  sections,  and  there  are  designs 
under  construction  which  do  not  even  consider  a  passive 
survival  of  a  structure.   These  proposals  include  a  vibra- 
tory icebreaker  motion  [7]>   ice  melting  by  internal  heating 
[17] *  air  cushion  mobile  platform  [22],  construction  of 
protective  soil  berms  [9]>  or  even  mobile  platforms  that 
only  drill  in  optimum  conditions. 

Scandinavian,  Soviet,  and  Argentinian  oil  interests 
are  not  represented  in  the  literature. 
1.5  Historical  Treatment 

Analytical  treatment  of  ice  pressures  and  winter  loading 
was  first  presented  in  fresh  water  civil  works,  especially 
dams  and  bridge  piers.   Early  investigators  faced  the  problems 
by  correctly  deducing  that  no  matter  what  the  loading  condi- 
tion or  geometric  possibilities,  ice  strength  could  not  exceed 
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its  crushing  strength.   Thus  a  limiting  value  of  Ocr  was 
imposed  into  the  following  relationship: 

FMAX  =bh  6CR  (Li) 

b  represents  pier  width  and  Ocr  crush  strength  of  ice 
'  Max  -*-s  force  exerted  on  structure 

A  special  study  of  the  A.S.C.E.  in  1931  (Committee  on 
Power  Division)  recommended  a  crush  strength  of  400  psi  [32] • 
A  review  of  Russian,  Polish,  Canadian,  and  U.S.  military 
design  codes  was  conducted  by  Watts  [32],  and  it  reveals 
that  most  design  criteria  are  simply  a  refined  empirical 
version  of  (1.1). 

Design  criteria  for  offshore  structures  in  the  much 
more  quantitatively  unknown  sea  environment  do  not  simplis- 
tically  prescribe  this  relation.   The  paragraph  relating  to 
design  for  sea  ice  forces  for  fixed  offshore  structures  from 
Pet  Norske  Veritas  1974  Rules  states  [10] : 

B401.   For  structures  intended  to  be  installed  in  areas 
where  ice  hazards  may  exist,  relevant  statistical  data 
for  the  area  in  question  are  to  be  submitted.   The  ice 
conditions  are  to  be  described  with  particular  attention 
to: 

-  concentration  and  distribution  of  ice 

-  types  of  ice  (ice  floes,  ice  ridges,  rafted  ice,  etc.) 

-  mechanical  properties  of  ice 
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-  mean  thickness  of  ice  floes 

-  average  drifting  speed  and  direction  of  ice 

-  probability  of  encountering  icebergs 

-  tidal  range 

This  rather  nebulous  requirement  is  in  stark  contrast 
to  the  curves,  tables,  formulae,  and  data  presented  for  other 
force  loadings.   It  is  here  that  the  lack  of  analysis  in 
the  state  of  the  art  is  most  keenly  felt. 

To  increase  the  body  of  knowledge  presently  available 
in  the  science  regarding  ice  pressure  loadings,  numerous 
investigators  in  the  last  decade  have  proposed  many  alternate 
procedures.  The  more  promising  approaches  will  be  analyzed 
in  a  subsequent  chapter. 

The  proposed  solutions  to  the  ice  pressure  problem 
fall  basically  into  four  broad  categories: 

1)  Classical  elastic  thin  plate  theory 

2)  Finite  element  or  finite  difference  numerical  solutions 

3)  Empirical  formulation 
K)     Model  basin  testing 

Each  approach  has  merits  as  well  as  severe,  restrictive 
limitations.  For  instance,  the  classical  theory  assumes 
thin  plate  elastic  behavior  (plane  stress  or  plane  strain, 
isotropic,  homogeneous,  elastic)  which  is  most  certainly  not 
the  complete  behavior  of  sea  ice.   Yet  the  solution  can  yield 
a  fair  degree  of  accuracy  within  certain  limiting  conditions 
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and  provide  a  closed  form  expression  which  is  easy  to  use 

and  easily  reproducible. 

1.6  Importance  of  Ice  Forces 

To  insure  that  the  reader  is  qualitatively  aware  of  the 
enormous  magnitude  that  ice  pressures  can  exert,  a  simple 
model  is  presented  in  Figure  1-1. 

A  monopod  platform  with  a  radius  at  ice  of  8  feet  is 
subjected  to  an  ice  load  driven  by  wind  or  current  shear.   If 
the  extent  of  the  floe  is  great  enough,  it  is  obvious  that 
forces  well  in  excess  of  local  ice  failure  can  be  generated. 
Using  a  (bCfl  of  200  psi  and  formula  (1.1),  adjusted  for 
projected  frontal  area,  the  ice  force  can  be: 


FMAX^  6cRbh  -  690  tons 


MAX     ^CR 

This  approximately  700  tons  is  concentrated  on  a  three 
foot  slice  of  the  leg,  presenting  severe  structural  consid- 
erations.  If  this  structure  were  located  in  80  feet  of 
water,  the  overturning  moment  could  be  20,700  foot -tons,  a 
serious  concern  to  foundation  and  soil  resistance  as  well  as 
flexural  strength  of  the  leg. 

Perhaps  an  even  more  serious  consideration  which 
has  recently  been  addressed  [21  ]  [25]  [26]  is  the  fact 
that  this  force  can  thrust  at  a  frequency  close  to  one  hertz. 
Apart  from  resonant  considerations,  the  fatigue  strengths 
of  the  cold  metal  and  cold  weldments  become  of  serious 
concern. 
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FIGURE    1-1 


Ice  interaction  with  monopod  type  structure 
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Peyton  [26]  has  summarized  some  important  requirements 
for  the  ice  force  design  of  marine  structures.   These  include 

-  Horizontal  forces  due  to  crushing,  tending,  buckling, 
or  shearing  of  the  ice  sheet  (quasi-static) 

-  Impact  forces  (dynamic) 

-  Weight  of  attached  sea  ice  at  lo\i   tide 

-  Buoyant  lift  of  above  at  high  tide 

-  Vertical  components  of  sheet  failure 

-  Diaphragm  bending  during  water  level  change 

-  Ice  accretion  by  spray 

-  Thermal  expansion  of  ice  in  joints 

-  Rubble  in  structural  framing  carried  by  ice 

-  Abrasion  by  moving  ice 

It  can  be  seen  that  in  a  cold  environment,  ice  forces 
can  be  large  and  obvious  as  well  as  an  insidious  problem. 
In  the  next  section  the  enormity  of  the  design  equation 
and  the  limits  of  present  knowledge  are  presented. 
1.7  Ice  and  the  Enormity  of  the  Design  Problem 

Ice  is  a  naturally  occurring  material  which  forms  in  a 
convenient  plate  shape  located  at  the  surface  of  the  liquid- 
air  interface.   Ice  formed  on  a  small  fresh  water  pond  can 
be  clear,  homogeneous,  constant  in  thickness  and  have  a  sen- 
sible temperature  distribution.   In  this  case  the  engineer 
could  possibly  apply  some  of  the  more  well-known  classical 
plate  solutions  to  find  the  ice  behavior  and  predict  with  a 
fair  degree  of  confidence  any  desired  solution  to  his  prob- 
lem. 
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Sea  ice  is  not  that  well  "behaved,  however.   Consider 
the  situation  shown  in  Figure  1-2  which  depicts  a  profile 
of  arctic  coastal  ice.  A  bottom  founded  structure  could 
be  located  anywhere  in  the  picture. 

Figure  1-3  presents  a  careful  summary  of  the  unknown 
parameters  affecting  the  design  equation.  A  high  degree 
of  interdependence  is  characteristic. 

Figure  1-4  presents  four  logical  paths  to  the  final 
design.  All  have  been  used  in  successful  designs.  Again, 
each  path  has  its  advantages  and  disadvantages. 

1)  From  known  environmental  conditions  -  this  is 
perhaps  the  ideal  situation;  a  successful  design 
could  possibly  he  created  anywhere  at  minimal 
cost,  once  perfected.   Of  course,  the  lack  of 
knowledge  regarding  the  environment  at  each  locale, 
the  rudimentary  state  of  analytical  tools,  and  the 
lack  of  correlation  with  known  quantities  makes 
this  procedure  years  away,  and  perhaps  unattainable. 

2)  In  situ  -  This  will  provide  better  data,  hut  still 
suffers  from  analytical  model  imperfection,  as  above. 
In  addition,  the  gathering  of  this  data  can  he 
costly. 

3)  Near  full  scale  -  a  very  practical  approach,  but 
it  too  can  be  costly,  nearly  as  much  as  the  final 
product  in  some  cases  [25].   It  does  offer  the 
advantage  of  risking  a  small,  inconsequential  failure 
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Parameters  affecting  structural  load  in  marine  conditions 
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over  a  larger,  more  serious  accident. 
k)     Previous  experience  -  This  approach  would  yield 
the  highest  confidence  levels  but,  unfortunately, 
little  has  been  done  and  some  of  what  has  is 
proprietary  information.  Here  again  an  overly 
conservative  design  would  propagate  other  overly 
conservative  designs. 
Without  presenting  a  complete  review  of  the  parameters 
in  Figure  1-3  or  the  literature  on  Figure  1-4,  both  of  which 
are  beyond  the  scope  of  this  thesis,  it  is  hoped  that  an 
appreciation  of  the  enormity  of  the  design  effort  involving 
ice  can  be  gained. 
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CHAPTER  2 
Summary  of  Some  Sea  Ice  Material  Properties 

2.1  Purpose  of  this  Chapter 

The  purpose  of  this  chapter  is  to  present  a  brief 
outline  of  some  important  sea  ice  parameters  and  properties 
for  use  in  developing  the  analytical  model  presented  later 
in  this  thesis.  For  a  thorough  treatment  of  ice  properties, 
the  interested  reader  is  referred  to  references  [12],  [13], 
[18],  and  [33]. 

2.2  Problem  Formulation 

Sea  ice  can  be  considered  to  be  sensitive  to  temperature, 
salinity,  crystal  orientation,  and  seasonal  variation.   It 
can  be  considered  to  behave  elastically  for  a  short  load 
duration  and  a  moderate  load  magnitude.  Katona  and  Vaudrey 
[18]  describe  an  approach  to  an  appropriate  field  theory 
which  is  sufficiently  general  to  apply  to  the  case  of  a 
circular  pile  -  ice  sheet  interaction  problem.  The  basic 
parameters  are  presented  in  Figure  2-1. 

Katona  and  Vaudrey  further  reduce  the  list  of  material 
parameters  (category  I  of  Figure  2-1)  to  a  functional  form 
which  contains  only  two  parameters : 


t 

(2.1) 


<V  frjKi[eKI(t/),x],-t/ 

where  Q..     is  the  stress  tensor  at  spatial  point  X 
and  time  T 
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FIGURE  2-1 

parameters  to  define  an  appropriate  field  theory 
(after  Katona  and  Vaudrey  [l8j  ) 
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r; 


k)  is  a  symmetric  fourth  order  tensor  functional 

allowing  complete  anisotropic  formulation 

£  (■{;')  is  ^ae   strain  tensor  history   0<"t  <  "t 

X    is  the  spatial  point  included  for  any  variation 

of  other  parameters 

X    current  time 
By  exploiting  knowledge  of  sea  ice  behavior,  (2.1)  can  be 
reduced  to  a  viscoelastic  representation  by  restricting 
load  magnitude  and  load  duration.  Further  restriction 
results  in  an  elastic  formulation.  The  resultant  material 
behavior  can  be  characterized  as  falling  somewhere  into 
Figure  2-2.   The  quantitative  limits  of  the  regions  have 
not  yet  been  experimentally  determined. 
2.3  Failure 

Category  III  of  figure  2-1  requires  knowledge  of  the 
mode  of  failure  of  ice  surrounding  a  circular  pile.  Most 
investigators  work  with  an  assumption  of  compression  failure 
of  the  ice  immediately  preceding  the  pile.  However,  since 
this  failure  is  quite  large  with  respect  to  shear  or  tensile 
failure  in  simple  tests,  it  is  conceivable  that  tension  or 
shear  failure  may  precede  compressive  failure. 

Additionally,  if  it  is  assumed  that  the  greatest  level 
of  force  on  the  pile  will  be  produced  at  the  point  of  rupture, 
cracking  and  large  deformation  behavior  need  only  be  applied 
for  postfailure  behavior. 

Failure  by  limiting  strain  is  possible  also  at  levels  of 
force  below  elastic  limits  and  at  larger  time  durations. 
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FIGURE  2-2 

Effect  on  domain  boundaries  by  load  parameters  at  a 

constant  temperature  (after  Katona  and  Vaudrey  P)8  ) 
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2.4  Material  Properties 

There  is  a  large  gap  in  the  body  of  knowledge  regarding 
the  material  properties  of  ice  and  nowhere  is  there  a  con- 
sistent and  universally  recognized  presentation  of  values. 
Therefore,  for  purposes  of  construction  of  the  model  and 
presentation  of  results  of  this  thesis,  non-dimensional 
quantities  will  be  used  wherever  possible.   Some  representative 
quantities  which  may  be  applied  for  clarity: 

Quantity  Typical  Value 

Young's  modulus  ErSp      300,000  psi 

Pois son's  ratio  H        .33333 

Shear  modulus  577+vO    112,500  psi 

Viscous  constant  «=*  ref      1.2788  x  10^  (psi-sec)" 

Elastic  limit  C>^p      150  psi 
Failure  by  simple: 

•compression  400  psi 

•shear  117  psi 

•tension  200  psi 

These  values  are  not  considered  to  be  rate  or  time 
dependent  for  this  model  formulation,  although  in  general 
they  are. 
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CHAPTER  3 
Analytical  Modelling 

In  the  existing  literature,  two  basic  mathematical 
models  were  found  for  the  ice  sheet  -  circular  pile  prob- 
lem. Both  models  will  be  quickly  reviewed  here  and  the 
features  and  limitations  discussed. 
3.1  Model  of  Frederking  and  Gold 

This  model,  presented  by  reference  [11],  draws  its 
basis  from  the  earlier  mathematical  treatise  of  Noble  and 
Hussein  [24].   In  this  earlier  work,  Noble  and  Hussein 
derived  the  exact  solution  for  dual  trigonometric  (Fourier) 
series  that  arise  in  the  solution  of  a  mixed  boundary  value 
problem,  such  as  that  one  depicted  in  Figure  3-1. 

The  assumptions  and  limitations  inherent  in  the  analysis 
of  Frederking  and  Gold  are : 

-  Interface  friction  is  not  considered 

-  Plain  strain  is  assumed 

-  An  infinite  elastic  medium  is  used 
Unfortunately,  the  series  presented  as  the  solution  is  red- 
ucible only  for  the  special  case: 

-^(l-P1)  =  (I-**)  (3.1) 

(primed  quantities  refer  to  properties  of  the  pile) 
Frederking  and  Gold,  however,  departed  from  this  analysis  early 
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FIGURE    3-1 

Indentation  geometry  (after  Prederking  and  Gold  \\\\   ) 
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and  used  the   stress   function     $>    suitable  for  the  boundary 

<t>     -,TT(i-p)L2(l-^)v©sme  -  (UZv»)[rUv-  + 
"£7]cG$QJ+D0Ur  ^|7*"*^   "**  "^J  D^e^G    (3.2) 


conditions j  r 

~  r: 


r1" 


h«2 

[where  r^  is  total  force  on  the  pile  concentrated  at 

the  center  of  the  rigid  pile] 
to  derive  the  radial  stress  and  strain.   The  infinite  series 
of  (3.2)  was  then  greatly  simplified  by  assuming  a  Poisson 
ratio  of  .5*  whereupon  the  representation  for  radial  strain 
at  the  interface  becomes: 

2Gf,Ke)-  -Sl5SES  +  £l  (3.3) 

The  constant  D0  can  be  uniquely  determined  in  terms  of  F^  . 
Frederking  and  Gold  then  go  on  to  develop  a  refinement 
of  (3.3)  in  which  temperature  effects,  viscoelasticity,  and 
a  relation  between  radial  strain  and  penetration  rate  are 
discussed.  Finally  an  equation  relating  total  force: 


Fy  =  K  (T)  L  (©)  6, 


t* 


clVi  (3.4) 


where  K(T)  is  a  temperature  correction  function 

$   is  a  constant  modifying  yield  stress  by 
strain  rate  effect 
and  L(6)is  a  geometry  function  based  on  (3.3)  given 
by: 
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LCe)  --      (cosg)  de  (3.5) 

here  ?}     is  the  "contact"  angle 
3.2  Model  of  Ross,  Hanagud,  and  Sidhu 

Ross,  et  al.,  [27]  considered  adaptation  of  one  of  their 
earlier  works  [34]  concerning  elastic-plastic  plate  impact 
study  to  the  problem  of  an  ice  floe  surrounding  a  rigid 
inclusion.  The  model  chosen  for  this  study  is  based  on  a 
finite  difference  scheme  for  the  solution  of  a  plane-stress, 
dynamic  problem.  The  problem  is  formulated  with  a  two 
dimensional  cartesian  network. 

The  method  presented  contains  the  following  assumptions, 
although  it  is  clear  that  the  model  is  not  restricted  to 
these: 

-  pile  is  rigid 

-  loading  is  by  edge  dislocation  or  stress  on  a  finite 
plate 

-  perfect  adhesion  of  ice  to  the  pile 

-  plane  stress  situation 

The  method  of  solution  of  this  dynamic  problem  is 
presented  in  reference  P4],   It  consists  basically  of  solving 
the  continuity  equations  for  the  body  forces  and  hence  particle 
accelerations.  The  accelerations  are  translated  into  strain 
rates  through  kinematic  relations.  Finally  the  stress -strain 
relations  are  used  to  determine  stress  rates.   These  stress 
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rates  are  then  integrated  numerically  into  new  stress  quantities 
whereupon  the  process  is  repeated. 
3.3  Criteria  for  the  new  model 

For  the  developement  of  a  model  for  analysis  of  the 
circular  pile  problem,  it  is  felt  that  the  following  items 
should  he  incorporated  into  the  ideal: 

-  The  analysis  diould  present  variables  to  include 
stress,  strain,  their  rates  and  displacements 

for  the  complete  spatial  distribution  and  at  various 
time  intervals. 

-  The  model  should  at  least  be  a  two  dimensional 
representation  and  be  capable  of  plane  stress  or 
plane  strain  variation. 

-  The  model  should  be  flexible  enough  to  accept  the 
various  material  representations  of  ice  including 
non-linear  and  time  dependent  ones. 

-  Geometrically,  the  model  should  be  capable  of  handling 
mixed  boundary  values,  interface  friction,  and  finite 
or  (semi)  infinite  sheet  formulation. 

-  Fracture  or  material  failure  can  be  easily  be  incor- 
porated. 

-  The  model  must  of  course  be  simple  to  use  and  econ- 
omical to  apply. 

Unfortunately,  these  ideals  are  often  mutually  exclusive. 
In  the  next  section  and  in  the  subsequent  chapters,  a  model 
will  be  developed  that  meets  at  least  several  of  the  criteria 
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and  will  "be  shown  to  be  useful  for  application  to  several 
ice  -  circular  pile  interaction  problems. 
3.4  Developement  of  the  New  Model 

The  model  to  be  developed  in  chapter  four  is  a  departure 
from  the  one  described  by  Ross,  et_al . ,  in  section  3-2.   The 
approach  used  by  Ross  seemed  interesting  in  that  it  is  capable 
of  meeting  several  of  the  criteria  required. 

The  method  of  solution  remained  critically  dependent  on 
the  time  integration  procedure,  however,  and  any  numerical 
inaccuracies  or  instabilities  would  propagate  throughout  the 
solution.  The  dynamic  problem  also  seemed  inappropriate  for 
study  of  slowly  varying  stress  fields.   Lastly,  Ross  gives  no 
results  for  his  dynamic  problem. 

The  model  proposed  in  this  thesis,  although  similar  in 
conception  to  Ross*,  exploits  the  radial  symmetry  associated 
with  the  physical  problem  to  work  in  only  one  dimension. 
The  method  of  solution  is  greatly  dissimilar  in  that  it  does 
not  treat  the  dynamic  problem  and  uses  a  linear  system  of 
equations  to  arrive  at  the  final  result.  The  new  model 
also  treats  the  case  of  linear  viscous  behavior  and  has 
capabilities  of  extension  into  other  material  behavior. 
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CHAPTER  4 
Developement  of  a  Finite  Difference 
Elastic  Constitutive  Model 

4.1  General 

The  equilibrium  equations  "by  Wang,  reference  [31  ], 
are  applied  to  the  element  depicted  in  Figures  4-1  through 
4-4.   In  two  dimensional  polar  co-ordinates: 

a»-  +  T  Je~  +   — r —   ~  ^  C4-1) 

7  36  +  "dF"  +   7  bre     =  ^Ue  (4.2) 

The  right  hand  sides  are  body  forces;  the  dots  indicate 
differentiation  with  respect  to  time. 

Since  sea  ice  is  a  transversely  orthotropic  material, 
variation  of  mechanical  properties  in  the  Z  Co-ordinate 
direction  will  he  due  mainly  to  salinity  variation,  temper- 
ature profile,  and  non-uniform  thickness.   For  the  purpose 
of  this  study,  however,  the  ice  sheet  will  he  considered 
to  he  uniform  in  thickness  with  zero  temperature  gradient. 
(Extensions  of  the  present  model  to  account  for  these  effects 
are  possible.)  For  this  reason,  developement  of  (4.1)  and 
(4.2)  into  a  two  dimensional  plane  stress  model  is  an 
appropriate  choice. 

Before  proceding,  it  is  perhaps  worthwhile  to  state 
the  assumptions  of  a  plane  stress  formulation: 
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FIGURE    4-1 


Stresses  acting  on  a  small  element  ABCD  in  polar 
coordinates 
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FIGURE    4-2 


Variation  of  direct  and  shear  stress  with  angular  change 
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FIGURE  4-3 
Compatibility  relationship 


FIGURE  4-4 
Sign  convention 
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<5£)6r2)0ez  are  zero 'on  "both  faces  of  the  sheet 

-  Therefore,  O^  C^C^ are  zero  "throughout  the  sheet. 

-  The  state  of  stress  can  be  specified  "by  6r,6e,£>re 
only  and  variations  of  these  quantities  are 
independent  of  vertical  position. 

4.2  Developement  of  the  Plane  Stress  Elastic  Model 

After  Wang  [31],  the  strain-displacement  relations 
for  small  deformations  are: 


€w-i^ur+^Lu..JSf  (4.5) 


Stress-strain  relationships  dictated  by  the  plane  stress 
assumption: 
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Combining  (4.1)  through  (4.8),  the  following  five  equations 
are  obtained. 

1.  fi  +  _L  A_  ft         ■   6r  -  6e    n  (Z|.9) 


H«.46»4»r.  -O  (*-10) 


£"r  -y  (0,.-  H0e)  -O  t4-11) 

Fu>"  +-f-^Ue--i-(6e-  K6r)  -  0      (*.12) 
F^Wr+^-Ue-fa,-2^  6rs.O  (4.13) 

The  body  forces  expressed  in  (4.1)  and  (4.2)  were 
set  to  zero  since  it  is  not  the  purpose  of  this  study  to 
expose  dynamic  considerations  or  impact  shock.   For  the 
range  of  validity  of  this  model,  inertial  effects  are 
negligible. 

Equations  (4.9)  through  (4.13)  (or  their  cartesian 

counterparts)  could  be  attacked  directly  in  finite  difference 

form,  but  the  ensuing  number  of  equations  would  be  quite 

large  and  difficult  to  apply. 

If  boundary  conditions  are  not  to  be  mixed  at  a 
1 
continuous  boundary,  a  Fourier  simplification  is  possible: 


Let   6^,69  and  Uj-  be  expressed  as  a  Fourier  cosine 

(even)  series, 

Let   6y-e,He  be  expressed  as  a  Fourier  sine  (odd) 

series, 


4o 


CO 


Orfr,©)  *  ^>    '  6r>n(ir)  tos  mG  (4.14) 

wi  =  o 

G^,(jCr    are  similar 
oo 

6re(r>e)  =  ^>    ,'   6re>n(i0  sin™8  (4.15) 

Ue     is   similar 


It  is  then  possible  to  express  (4.9): 


00  oo 


M  2  6rTOWcosme]+  f^[£0UBWsmme]  + 
oo  oo 

w»=o  -»      '    L  '      J       m  J 

(4.16) 


Vw'O 


which  reduces  to: 


CO  co 


S^6^ cosine  +  f-]T]  ^cosTnG  + 

CO 

+  T  2i-j  (  °U  "    6<^  )  co*  ™6  (4.17) 

>ViaO 


Finally: 


^6r0+0  +  f(6re-66o)j  ,  o 
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(4.18) 


Thus  Of    can  "be  expressed  as  an  equation  relating  only 
like  wavenumbers.  The  set  of  solutions  unique  to  each 
wavenumber  can  then  be  reconstituted  into  a  Fourier  series 
by  using  ah  appropriate  number  of  waves. 

Equations  (4.10)  through  (4.13)  can  similarly  be  trans- 
formed into  an  infinite  number  of  sets  of  five  equations: 


Tf  °r*i  +  T  $™ 


yn 


4-  ^r-wi  -  6ewi 
in       f 


-T^+^Ove^r  ° 
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TOO 


^Ur*   -^(C)rw,-  W6ewi) 


w\   ^r  ^™   If   ew     r    ^r© 


r  Uirv  +^.  ue™  ~T"  U 


(^.19) 


=  0 


™J(O<Y*<00) 
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°©^ 


in  the  fourth  row  of  (4.19)  can  be  explicitly 
expressed  as: 


<V  -  |  urM  +  £*  Uew  +  y  o, 


W\ 


(4.20) 


Since  no  equation  contains   a  derivative   of     6q.       ,   this 
variable  can  be  eliminated,  yielding: 


Yf\ 
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■Ewi 


2>v- 


VC-^XV*  I'-f"  )6re+(^r)u, +(^f)ue 


£6~+F^)<V  +    (4-   )6re+(--^)U,+r-^-)u( 


Sr  MV  v       V 
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WW 


£uf+(-^-)<5r  +  {   o    )6re+(-^-)Ul.+(-i:r-)ae 


We     O      )0r+    r^JCrt  +  l^Ur^l-T-JU, 


Let  the  vector     r    represent: 


(4.21) 


=  O 
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lA4 


Oim  £  co 


(4.22) 


Then  (4.21)  becomes  the  set  of  linear  equations: 


^MZ-[A]Ztt=Q     <*•«> 


where  Fj J  is  the  4x4  identity  matrix  and 


where 


[A] 
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is  expressed  by  the  equation; 


[A]    = 


\-v> 
Y 

r 
•(1-^) 

0 


t 

2_ 

Y 


-E 


m 


0 

•2(I+h) 


YY\£ 


V_ 
r 

-m 
r 


-Em2 


(4.24) 


(4.23)  represents  an  infinite  set  of  4  simultaneous  first 
order  linear  differential  equations.   It  can  he  seen  that 
two  spatial  variables  would  result  in  an  equation  relating 
a  derivative  in  0  as  well. 

The  system  of  equations  expressed  by  (4.23)  can  only 
be  solved  numerically. 
4.3  Numerical  Solution  of  the  Plane  Stress  Elastic  Model 

The  method  of  finite  differences  will  be  exploited  to 
solve  equation  (4.23) 

Considering  the  finite  mesh  scheme  presented  in  Figure 
4-5  and  in  Figure  4-6,  the  appropriate  difference  equations 
valid  at  the  nodes  are  expressed: 

^i](-3Z,+  4Z„-y4A]Z„=Q 

valid  at  node  1  (4.25) 
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FIGURE  4-5 
Finite  difference  representation 


(i  o o o  o e oooo a  t 


12       3       4 


1-2    1-1    I 


FIGURE  4-6 
Node  numbering  scheme 
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valid  at  the  I-th  or  last  node, 
and  the  central  difference 


3fc[l]f  Z.-ZJ+[A]Z.  =  0       (*.*> 


valid  at  all  interior  nodes. 
Recasting  equations  (4.25),  (4.27),  and  (4.26): 


1   {M-m]±,+mz«mz,--Q 

•   ,         •  •  (^.28) 

^WZ.-KilZ^ftAl^LiljZ/O 

(4.30) 

Equations  (4.28),  (4.29),  and  (4.30)  are  then  loaded 
into  a  suitable  simultaneous  equation  -  ,matrix  representation 
as  shown  in  Figure  4-7. 

This  representation  is  then  expanded  into  the  full  coeff- 
icient matrix  shown  in  Figure  4-8.   It  should  he  noted  that 
even  the  simplest  scheme,  a  three  node  model,  will  require 
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reduction  of  a  144  element  square  matrix  (12  equations  in 
12  unknowns).  Even  a  numerical  scheme  such  as  the  Gauss 
elimination  or  the  Gauss-Jordan  reduction  scheme  is  too 
laborious  for  hand  calculation,  so  a  machine  solution  must 
be  employed. 

Figure  4-8  represents  the  expansion  of  equations  (4.28) 
through  (4.30).  The  elements  noted  represent  the  corres- 
ponding position  in  (4.24).  As  an  example: 


A„=f 


r  is  evaluated  at  node  2 


Elements  labeled  R  or  £..  represent  modification  to   Al 
required  by  forward  difference  and  backward  difference 
equations  (4.25)  and  (4.26).  The  matrix  is  a  square  matrix 
distinguished  by  a  band  twelve  elements  wide.  Even  for  a 
small  number  of  node  equations  the  ensuing  sparse  matrix  is 
wasteful  in  space  and  computation  time.  Consequently  the 
representation  in  Figure  4-8  is  further  compressed  to  a 
band  storage  array  depicted  in  Figure  4-9. 
4.4  Non-dimensionalization 

In  order  to  produce  array  elements  of  approximately 
the  same  order  of  magnitude,  (4.23)  was  converted  to  non- 
dimensional  form  by  the  following  conversion: 

on  stress   ■  "   =  §     (non-dimensional) 
Oref 

where    O    =  F"  „— 
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FIGURE  4-9 
Matrix  in  the  band  representation 
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is  taken  as  the  average  Young's  modulus  through 


the  thickness,  and  length  is  non-dimensional  by 

L 


L 


=   L 


REP 

where 


L  =  R 


IN 


(non-dimensional) 


or  inner  radius 


Equation  (4.23)  becomes  non-dimensional  and  the  coefficient 
matrix  (4.24)  becomes: 


Here 


l-v> 

Vr\ 

-1 

-Wl 

Rao 

RAO 

RAO* 

RAO2 

-Vr\lr> 

2 

-Wl 

2 

RAO 

RAO 

RAO1 

RAO* 

0 

RAO 

RAO 

0 

-2(l+v) 

RAO 

-I 
RAD 

-        1    +    /   L 

-ivRcwt-n 

(4.31) 


J-  (4.32) 

and  will  vary  within  the  array  according  to  node  position  i. 


4.5  Machine  Solution 

Equation  (4.23)  was  implemented  in  a  Fortran  language 
subroutine  called  elastic.  For  reasons  discussed  in  the  next 
section,  the  number  of  nodes  used  was  400.  Consequently  the 
core  storage  requirements  in  double  precision  amounts  to 
approximately  60K  words. 
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Solution  of  (4.23)  was  by  means  of  an  I.M.S.L. 
subroutine  called  LEQT1B  [16].   This  method  of  solution 
of  a  system  of  linear  equation  is  described  by  [23]  as 
the  Crout  factorization.  The  Crout  factorization  is  similar 
to  Gauss  reduction  methods  but  the  lower  triangular  decom- 
position is  used  instead.  Partial  pivoting  and  row  normal- 
ization are  employed. 

LEQT1B  was  slightly  modified  into  a  double  precision 
routine  and  labeled  dpband  by  the  author.  This  routine  is 
included  as  an  integral  part  of  subroutine  elastic  and  is 
listed  in  Appendix  A. 

Subroutine  elastic  was  written  for  the  multics  system 
implimented  on  a  Honeywell  6l80  used  by  the  Massachusetts 
Institute  of  Technology  Information  Processing  Center. 
Approximate  CPU  time  for  one  solution  was  16  seconds. 

A  flow  chart  for  subroutine  elastic  is  presented  in 
Figure  4-10.  A  symbol  table,  description  of  common  blocks, 
listing,  and  a  user's  guide  are  included  in  Appendices  A-l 
through  A -4. 
4-6  Numerical  Stability 

The  number  of  stations  chosen  is  400.  Table  4-1  and 
Figure  4-11  illustrate  the  convergence  of  the  numerical 
solution  as  the  number  of  stations  is  increased  for  the 
particular  case  of  zero  wavenumber.  Having  checked  the 
convergence  of  the  numerical  solution,  the  question  now  is 
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FIGURE  4-10 
Block  diagram  for  subroutine  Elastic 
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Table  4-1 


tfavenumber 


Maximum  Oscillation 
(per  cent) 


0 
1 
2 
3 
4 
5 
6 

7 

8 

9 
10 


0.0075 
1.7500 

0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 


Radius  -  outer  =20. 

Radius  -  inner  =  1. 

u  =  u„  =  0.   inner  boundary 

6r  =  -.66666667  outer  boundary 


6re=  0. 


outer  boundary 


Maximum  oscillation  in  numerical  solution  by  wavenumber 
(solution  carried  to  five  digits) 
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whether  the  numerical  scheme  is  consistent,  i.e.  whether 
this  convergent  numerical  solution  is  the  true  solution. 
This  is  done  in  the  following  section. 
4.7  Comparison  with  Analytic  Solution  and  other  Tests 

Fortunately,  the  analytic  solutions  to  the  relations 
expressed  in  equations  (4.21)  can  be  obtained  and  are 
presented  in  Appendix  B.  The  analytic  solutions  and  their 
computed  counterparts  are  illustrated  in  Figures  4-12 
through  4-14  for  wave  numbers  0,  1,  and  2  respectively 
for  the  case  of  radial  stress  at  R  out/R  in  =  20. 

Test  A  This  was  run  to  determine  the  effect  of  radius 
on  the  inner  stresses.   It  would  be  expected  that  for  a 
stress  loading  of  the  outer  boundary,  the  importance  of 
the  radius  on  the  inner  boundary  would  decrease  with 
increased  plate  diameter.  This  has  been  determined  to  be 
so  numerically,  and  some  results  are  presented  in  Figures 
4-15,  4-16,  4-17,  and  4-18.   It  is  also  noteworthy  that  the 
higher  wave  numbers  have  very  little  effect  on  the  inner 
stress  solution.  A  conclusion  here  can  be  drawn  that  will 
save  much  computational  effort  if  only  a  limited  amount  of 
data  is  desired. 

Test  B  This  was  run  to  determine  the  sensitivity 
of  the  solution  to  the  Poisson  effect.   The  results  for 
four  wave  numbers  are  presented  as  Figure  4-19.  Wave 
number  one  appeared  to  be  almost  entirely  insensitive, 
whereas  wave  number  zero  appeared  to  be  the  most  sensitive. 
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6r=-666667 
ROUT =20 


WAVENUMBER 
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POISSON  RATIO '^ 


Variation  of  radial  stress  at  the  inner  boundary 

as  a  function  of  poisson  ratio  for  several  wavenumbers 
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The  next  two  tests  were  of  load  magnitude  and  shape 
and  designed  to  test  for  linearity. 

Test  C  The  load  magnitude  was  varied  by  integral 
i  multiples,  retaining  the  same  distribution.   Results  were 
exact  multiples  as  well.  No  results  are  presented. 

Test  D   The  test  of  load  shape  was  more  of  a  test 
of  the  Fourier  subroutine  which  converts  boundary  condi- 
tions into  series.   The  following  radial  loading  was 
supplied : 

&  #2  #3 

Or     a  1.cos6  a  1.-1.  cose  «1. 

Ore  =  0.  =0.  -  =0. 

Ur  i«h«  =  0.  =0.  f  =0. 

U©  IMMSU  =  O.  =0.         =       -0. 

The  results  were  excellent  and  are  presented  in  Table  4- II. 

Test  E  The  last  test  concerned  the  radius  at  which 
the  finite  plate  could  be  considered  infinite  for  practical 
purposes.   This  test  is  similar  to  Test  A.   This  test  was 
somewhat  inconclusive  since  above  a  certain  radius  some 
solutions  became  unstable  again  and  innacurate.  A  practical 
limit  of  radial  spacing  was  chosen  at  n=400,  since  the 
marginal  accuracy  of  finer  differences  was  outweighed  by 
increases  in  computational  time,  cost,  and  round  off  error. 
Results  of  this  test  are  presented  as  Test  E  in  Table  4- III. 
4.8  Comparison  with  Ross'  Solution 

In  Section  3.2  a  review  of  the  analysis  of  Ross  was 
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Table  4-1 I 


Radius 


Run  1 


1.00 

9.72992 

2.90 

5.34873 

4.80 

3.36513 

6.70 

2.46372 

8.60 

1.95948 

10.50 

1.64134 

12.40 

1.42509 

14.30 

1.27047 

16.20 

1.15562 

18.10 

1.06816 

20.00 

1.00002 

urJue=0 

inner 

outer 
6r=COS0 

6©=0 
Rout  =20 


Run   2 

-8.2323 
-4.2905 

-2.3445 
-1.4538 

-0.95387 

-0.63805 

-0.42309 

-0.26917 

-O.I549I 

-0.067790 

-0. 11921 e-04 


urjue=0 

inner 

outer 
(br=1-cos9 

6e=0 

Rout=20 


Sum 

1.49762 
1.05823 
1.02063 
1.00992 
I.OO56I 
1.00329 
1.00200 
1.00130 
I.OOO65 

1.00037 
1.00000 


Run  3 

1.49759 

I.05824 
1.02047 
I.OO989 
I.OO551 
1.00328 
1.00200 
1.00119 
I.OOO65 
1.00027 
1.00000 


urju©  =  0 

inner 
outer 

6r=1 

6e=0 
Rout  =  20 


Test  of  elastic  routine  by  variation  of  load 
shape 
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Table  4-III 


tfavenumber 


3 
>4 


^T  outer 
^P  inner 
©Rout 


=10 


5%  deviation  Instability 
©Rout        ©Rout 


>100 

>100 

15 

>100 

50 

10 

>100 

>100 

75 

25 

>100 

>100 

<25 

>100 

>100 

Practical  extent  of  finite  condition  at  certain  wavenumbers 
by  variation  of  outer  radius 
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presented.   This  paper,  reference  [27] >   also  included  an 
example  in  which  the  loading  depicted  in  Figure  4-20  was 
analyzed.  Below  the  stated  elastic  limit,  it  was  possible 
\   to  compare  the  model  developed  in  this  chapter.  An  effort 
was  made  to  correlate  this  result  with  a  similar  loading 
applied  to  the  method  developed  in  this  chapter.   The 
comparison  is  presented  as  Figures  4-21,  4-22,  4-23. 

Although  the  basic  shape  and  critical  values  of  cross- 
over are  similar,  discrepancies  were  noted  with  respect 
to  the  hoop  stress  and  also  with  respect  to  the  order  of 
magnitude.  By  changing  the  displacement  loading  magnitude 
to  .01  inches  vice  .01  feet,  the  numbers  were  much  alike. 
No  explanation  was  found  for  the  difference  in  hoop  stress, 
but  from  equation  (4.20): 

Evaluated  at  Rw«  U©,Ur  =  0    (boundary  conditions) 
produces 

or  a  non-zero  value  for  any  non-zero  radial  stress. 

This  loading  was  also  run  with  3,  5*  and  11  wave  numbers 
and  no  effect  was  found  on  the  value  of  inner  face  stress, 
although  significant  variations  existed  toward  the  edge. 
This  is  as  expected,  since  the  influence  of  the  trigonometric 
components  vw)|  do  not  substantially  contribute  to  the  center 
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Free  Bounda::y 


12ft. 


Bouncary 


FIGURE    4-20 

Discontinuous  edge   loading  of  a  finite   plate   about 
a  rigid  inclusion     (after  Ross,   et  al,   [27] ) 
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stress  as  the  edge  of  the  plate  is  extended  toward  (perceived) 
infinity. 

Figure  4-24  displays  the  radial  loading  represented 
by  three  Fourier  coefficients  and  the  original  loading.   It 
can  be  seen  that  even  for  this  discontinuous  loading,  three 
wave  numbers  presents  a  good  approximation. 
4.9  Peculiarities  and  Range  of  Validity 

This  elastic  formulation  is  derived  from  the  basic 
plane  stress  equations  and  consequently  should  be  useful 
where  plane  stress  is  appropriate.  The  model  achieves 
significant  accuracy  up  to  a  ratio  of  outer  radius  to  inner 
radius  of  25:1  •  Beyond  this,  the  solution  diverges  from 
the  analytical  derivation. 

The  worst  case  is  at  wave  number  one  (Figure  4-25). 
It  is  doubtful  that  the  infinite  case  can  be  modeled  by  a 
finite  difference  technique  in  this  problem,   since  a 
boundary  loading  for  the  first  wave  number  at  all  radii 
will  produce  a  significant  solution  at  the  pile.   Consequently, 
this  elastic  model  is  useful  for  radial  limits  of  outer  radius/ 
inner  radius  ratios  less  than  25. 

The  buckling  limit  of  a  thin  plate  is  an  important 
consideration.  The  analytic  derivation  of  one  buckling  load 
is  presented  as  Appendix  C.  The  equations  of  plane  stress 
I  do  not  allow  for  deflections  in  the  third  dimension,  therefore 
they  will  not  predict  the  failure  by  this  mode.  External  means 
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FIGURE   4-24 

Displacement  loading  (radial)  applied  to  boundary 
of  Ross'  solution  -  comparison  with  first  three 
wavenumbers 
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must  be  applied  to  the  program  to  preclude  this  occurance. 

The  plate  is  assumed  to  be  attached  to  the  rigid  in- 
clusion (frozen-in  problem).  The  partial  movement  of  either 
boundary  cannot  be  supported  by  this  method,  although  it  is 
conceivable  that  a  more  general  two-dimensional  model  could. 

Finally,  the  model  can  support  an  application  of  forces 
or  displacement  at  either  boundary,  which  can  in  turn  be 
discontinuous.  The  degree  of  refinement  is  limited  by  the 
wavenumbers  considered.  The  accuracy  is  limited  by  the 
circumferential  spacing  chosen. 
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CHAPTER  5 
Viscous  Extension  of  the  Model 


|  5.1  General 

Considering  the  simple  rheological  model  in  Figure  5-1 > 
the  following  differential  relations  hold : 

•  •  • 
£r  s  £r  vtscoos  +  (^elastic  (5.1) 

•  •  • 

6©    —       6^  V1SCOOS        +        fe  ELASTIC  (5.2) 

•  *  • 

fre=    6r©v»scoos    +     £re  elastic  (5.3) 

where   f  vuicooa  *  Action  (  6 )        (5.4) 
The  Maxwell  model  is  chosen  for  its  simplicity  and 
because  it  is  possible  to  generalize  into  a  non-linear 
creep  law.  Additionally,  some  investigators  have  modeled 
sea  ice  behavior  this  way  experimentally,  and  physical 
values  are  available  [33]. 
5.2  Incorporation  of  the  Viscosity  Relations 

The  stress-strain  relations  of  (4.6),  (4.7) »   and  (4.8) 
can  be  expressed  using  (5.1)*  (5.2),  and  (5.3): 

;        (r   *  o<(0r- POe)  +  |r  (6,  -  v^6e)     (5-5) 

(e  -  <*  (  Oe-HCv)  +  -1  (6e-  w6r)    (5.6) 

|        („-  &Lii±22  Ore  (5.7) 
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M  =  H-N 


' 7    relaxati 


on  of 


stress 


//  /  p. d. -permanent 

' J   '  deformation 

//  r.s. -recovered 

'/  strain 


FIGURE  5-1 

Model  for  Maxwell  liquid  (after  [8]  ) 
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These  relations  can  be  incorporated  into  a  series  of  equa- 
tions as  in  section  4.2  yielding: 

i^[l]Z.*[A]Zj-[E]Z.    <».«> 

Where:  O  <  \m  £  CO 

A    is  the  same  matrix  as  expressed  by  (4.24) 


And 


[E]  = 


2         2 
-to*     -Ewc( 


-2         r2  X 
-fc  WvW    -fc  to  ol 


SotQ  +  »Q   Evtt*(U-»>) 


(5.9) 


2*  (l  +  w) 


5.3  Solution 

Numerical  solution  of  the  series  of  differential 
equations  (5.8)  is  carried  as  before,  only  the  unknown 

m 

variable,   7   ,  represents  rate  at  time  t    •  The  value 

used  on  the  right  hand  side  of  (5.8)  is  the  elastic 
solution  for  j?      .  Note  that  boundary  conditions  of  load 
rate  must  be  employed.  The  solution  yields: 


MZ 


m 


A   CoiaSTA»aT 


(5.10) 
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This  differential  equation  in  time  can  "be  solved  using  the 
Euler  or  "one  step"  method  to  arrive  at  a  solution  for  the 
next  time  step,  where  the  process  is  repeated. 

f  t        t  +  A-fc 

Non-dimensionalization  with  respect  to  another  variable 

is  required.  The  variable  chosed  was  oi  ,  the  viscosity 

-I 
coefficient,  with  units  of  |_  T  M  • 

If  ©Preference  is  chosen,  and  mass  is  preset  by  consideration 

of  stress  non-dimensionalization  (Section  4-.^),  then  derived 

time: 

~r       ! (5.12) 

Therefore  any  value  of  non-dimensional  time  is  related: 

Trcf 

As  an  example,  using  typical  values  for  sea  ice  [  ], 

Erep  =  20SS  p«sf 

*ref    =  1.2788*  I0"S  (psf-  sec)"1 

implies 

E=  i. 

ana  «*  =    I.  Tnus         -pR6F  _  2607  see 

or   one  unit  of  non-dimensional  time 
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A  flow  chart  of  the  numerical  computer  solution  and 
time  integration  is  presented  as  Figure  5-2. 

5.4  Convergence  in  Time 
The  Euler  method  chosen  for  time  integration  assures 

convergence  if  the  time  step  is  small  enough.  Since  it 
is  not  possible  to  analytically  determine  the  solution,  as 
was  done  for  the  elastic  case,  one  approach  is  to  run  a 
general  solution  for  increasingly  small  time  steps  until 
no  numerical  improvement  is  achieved.  Unfortunately  the 
solution  time  can  be  prohibitive.  The  results  of  two  time 
convergence  test  series,  for  wave  number  zero  and  one,  were 
found  to  converge  at  a  time  step  of  .001  "]~  s   beyond  which 
significant  improvement  was  not  achieved.  The  time  increment 
to  achieve  this  numerical  stability  comes  out  to  be  a  constant 
independent  of  stress,  for  the  upper  bound  prediction. 

5.5  Range  of  Validity 
The  viscoelastic  formulation  presented  in  this  chapter 

is  not  valid  where  the  elastic  theory  is  inappropriate  or 
where  the  numerical  limitations  of  subroutine  elastic  do 
not  apply.  The  range  of  validity  is  further  restricted 
on  how  well  and  on  which  range  of  parameters  ice  can  be 
modeled  as  a  linear  viscoelastic  material. 
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Subroutine  Visco 
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Read  time  step 


(l 


Read  stop  time 


Read  time  solutions  to 
be  displayed 


Save  elastic  solution 

Load  modified  solution 

vector 

Load  rate  boundary  conditions 

Solve  system  of  equations 

for  rates 

Solve  differential  equation 

Increment  time  step 


Loade 


Boundb 


Lpband 


display  interval  reached? 


.•fritr3| 


Strain 


time=stop  time? 


FIGURE  5-2 

Block  diagram  for  subroutine  Visco 
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CHAPTER  6 
Analysis  of  Ice  -  Circular  Pile  Behavior  with  the  Model 

The  previous  chapters  have  dealt  with  the  need  for  an 
analytical  tool  and  the  construction  of  a  mathematical 
model.  Without  exhaustively  analyzing  the  many  combinations 
of  loadings  and  important  parameters  that  affect  the  circular 
pile  problem,  some  interesting  results  will  now  be  developed. 
6.1  Material  Behavior 

Katona  and  Vaudrey  [18]  have  presented  a  logical  seq- 
uence of  obtaining  ice  material  behavior  experimentally  for 
a  particular  geometric  configuration.   The  results  of  these 
data  are  then  fitted  into  curves  which  can  then  yield  a  range 
of  validity  of  ice  material  law  as  discussed  in  Section  2.3 
and  with  Figure  2-2. 

Of  course,  since  this  model  formulated  in  this  thesis 
responds  only  viscoelastically,  the  transition  from  visco- 
elastic  to  viscoplastic  behavior  cannot  be  obtained.  The 
model  can,  however,  define  the  elastic-viscoelastic  transition. 

As  an  example,  consider  the  loading  situation  depicted 
in  Figure  6-1.  Physically,  this  could  represent  the  rigid, 
frozen-in  pile  loaded  by  a  large,  thermally  expanding  ice 
sheet.  Neglecting  any  temperature  effects  to  viscosity  and 
elastic  modulus,  the  loading  curves  of  Figure  6-2  are  con- 
structed.  Stress -strain  diagrams  similar  to  Figure  6-3  show 
the  viscoelastic  creep  behavior.   The  apparent  relaxion  of 
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Viscosity  coefficient    oi.    -    1 . 

Elastic  Modulus         E  =  1 . 

Non-dimensional  quantities, 

see  text  in  Sections  4-4  and  5—3 


-      -3 
br  =.222X10 

tore=0 


FIGURE    6-1 

Circular  ice  sheet  -  rigid  pile  geometry  and  loading 
configuration  used  in  viscoelastic  behavior  analysis 

of  Section  6.1. 
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Variation  of  strain  at  pile-ice 
interface 
u-0  Loading  condition  of  Figure  6-1  applies 


20  _ 


10  _ 


5  - 


2  - 


0 


Outer  boundary  conditions 
Radial  load  applied 


(=>  =  -.10x10 


-2 


=  -.  222x1 0"3 


=-.111x10 


-3 


— r~ 
.1 


1 — 
.2 


T — 

.4 


— l 
.5 


Load  duration  -  non-dimensional  time 


FIGURE  6-2 

Viscoelastic  behavior  -  Effective  strain  vs.  load  duration 
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6eff  =  (3(b^be-brbe+3bre)r 
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Non-dimensional   time 
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FIGURE  6-3 


Effective  Strain  xlO  J         J_ 

2   2       2  2 

eeff=(er+ee+6ree+€re) 


Viscoelastic  behavior  -  stress— strain  diagram  showing 
creep  for  loading  condition  of  Figure  6-1,  at  the  pile- 
ice  sheet  interface. 
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elastic  modulus  is  then  limited  to  a  value  of  10$  decrease. 
This  then  becomes  the  transition  boundary  between  visco- 
elastic  and  elastic  behavior  of  Figure  2-2. 

In  the  example 3  loadings  of  those  shown  in  Figure  6-2 
are  used  to  define  a  region  of  validity  in  Figure  6-4. 

Using  actual  values  such  as  those  found  in  Section 
2-4,  it  can  be  seen  that  a  time  of  approximately  10  minutes 
defines  the  region  at  low  stress  levels  (below  yield). 

The  effective  stress  and  strain  plotted  in  Figure  6-3 

are  defined :       2     i    2.    2  i 

6eff=T(0*-  +  °©  "  0v6e  +  3  6re)   (6.1) 

For  this  loading  condition,  it  can  be  seen  that  a 
value  for  modulus  of  elasticity  is  about  one  half  the  uni- 
axial value.   This  will  vary  as  well  along  each  point  in 
the  radius  of  the  ice  sheet.   Interestingly  enough,  the 
relaxation  to  90$  of  its  former  value  occurred  at  the  same 
time,  despite  the  variation.   Under  a  more  general  loading, 
not  confined  to  the  zero  wavenumber,  it  is  likely  that  the 
material  will  behave  not  in  the  same  region  of  Figure  2-2 
over  its  entire  geometry. 
6.2  Under  Rate  Loading 

The  example  in  the  previous  Section  was  subjected  to 
rate  conditions  on  the  loading.   This  was  done  to  test  the 
capacity  for  this  parameter  in  the  model.  The  load-time 
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FIGURE  6-4 

Viscoelastic  behavior  -  Elastic-Viscoelastic  transition 
from  load  curves  of  Figure  6-2.  By  definition,  relaxation 
of  elastic  modulus  is  limited  to  10$. 
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curves  for  the  loading  of  Figure  6-1  are  presented  in  Figure 
6-5  for  three  loading  rates.   The  greatest  rate  is  very 
close  to  the  maximum  of  5000  psi/min  suggested  by  Katona 
and  Vaudrey  to  avoid  dynamic  effects.   No  general  conclusions 
are  drawn  here  since  the  size  of  the  initial  load  and  the 
various  rates  have  little  significance  without  a  specific 
physical  problem  to  simulate. 
6.3  Analysis  of  a  Friction  Coefficient 

The  loading  situation  in  Figure  6-6   was  applied  to  the 
model  to  determine  the  ratio  of  shear  stress  to  (compressive) 
radial  stress  at  the  rigid  pile  interface.   This  could  be 
construed  to  be  an  indication  of  a  static  type  of  frictional 
coefficient,  but  since  the  model  also  assumes  capability 
of  the  interface  to  support  a  tension,  in  this  region  "fric- 
tion coefficient"  would  be  meaningless. 

The  analysis  proceeded  with  those  values  that  this  edge 
loading  produced.   Defining  the  coefficient  of  static  friction 
as : 


M*  ■ 


•re 


(6.3) 


it  can  be  seen  that  this  coefficient  will  vary  about  the 


interface  at  least  to  the  point  at  ®~~Z~    where  &r   becomes 
tensile  and  a  coefficient  of  friction  becomes  meaningless 
(Figure  6-7). 

At  a  value  of   0  *  .31  radian  (given  this  loading), 

j^s  *  .3 
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Effective 

Strain       xlO'  . 

€eff=(€r+€e+eree+£re) 


.3        .4        .5 

Load  duration  —  time 
T 


FIGURE  6-5 

Visco— elastic  behavior  —  Cgff  vs.  time  curve  under 
three  loading  rates,  at  the  pile  -  ice  sheet  interface. 
Loading  condition  of  Figure  6-1  applies,  with  6p=.llllxl0' 
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Rout  =20. 

2 
£)p  =-.333  cos  ©      outer  boundary 

fo~i=-.333  sin©  cos©  outer  boundary 

for  -tt  <  e  <  "]T 

and  0.  otherwise 


at  the  inner  boundary  or  pile-ice 
sheet  interface: 

Rin  =  1. 

Up  =Ufi  =  0.   considering  "friction" 

and 

bpQ  =  Up  =0.    "frictionless" 


FIGURE  6-6 

Loading  condition  circular  pile  -  ice  sheet  interaction 
for  cases  discussed  in  Sections  6-3  and  6-4. 
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Angular  position  in  radians 


frictional  considerations  valid 
in  this  region 


FIGURE  6-7 

Ratio  of  shear  stress  to  radial  stress  at  the  circular 

pile  -  ice  sheet  interface,  loading  of  Figure  6-6  applies. 
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Now  if  there  is  no  ice  -  pile  adhesion  present 
and  if  the  frictional  coefficient  of  smooth  sea  ice  on 
a  smooth  steel  circular  pile  is  about  .3  as  some  inves- 
tigators indicate,  then  it  is  clear  that  for  regions 
beyond  .31  radian  CC©  will  be  nonzero.  A  situation  in- 
volving mixed  boundaries  and  separation  will  than  be  evident. 

By  changing  the  loading  to: 

6>-  -   -  Oy  co*  0 
Or©  =  "  6r©  sm  G 

it  was  found  that  the  coefficient  of  friction  at  the  inter- 
face changed  only  slightly.   Since  for  the  wavenumber  zero, 
only  normal  stresses  are  conveyed  (being  an  odd  function) 
and  the  higher  wavenumbers  lose  their  influence  with  increased 
radius,  frictional  forces  are  important  only  in  the  region 
0|  £..31  radian  for  non-axisymmetric  loads. 

A  conclusion  can  be  drawn  that  for  the  model  of  the 
circular  pile  -  ice  sheet  interaction  and  unless  ice  adhesion 
is  present,  frictional  behavior  is  confined  to  a  small  area 
preceding  the  pile  and  may  be  unimportant  in  analysis. 
6.4  Analysis  of  the  Behavior  of  Ice  Floe  -  Pile  Interaction 
Comparing  Adhesive  and  Non-adhesive  Behavior 
Since  the  influence  of  friction  at  the  circular  pile  is 
minimal,  it  is  necessary  now  to  compare  the  behavior  at  the 
interface  using  complete  adhesion  (frozen-in  situation)  and 
and  using  complete  circumferential  freedom. 
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The  loading  of  Figure  6-6   again  was  applied  to  the 
model  to  consider  the  effect  on  the  radial  stresses  of  a 
"frictionless "  boundary.   Since  the  loading  is  the  same 
for  both,  the  net  force  felt  by  the  pile  is  the  same  in 
both  cases.   Figure  6-8  presents  the  variation  of  radial 
stresses  in  both  cases.   The  point  where  adhesion  is 
required  is  precessed  to  1.48  radians. 

Noble  and  Hussein  [24]  obtain  a  value  of  1.24  radians 
analytically  in  their  formulation.   Severe  differences  do 
exist  in  the  problem  as  their  model  does  not  support  tension. 
Therefore  it  would  be  expected  that  in  the  new  model, 
the  radial  stress  on  the  compression  side  would  be  supporting 
the  lost  force  on  the  tension  side,  with  basically  the  same 
shape  in  the  loading  curve. 

At  any  rate,  a  conclusion  can  be  drawn  that  the  friction- 
less  case  is  the  worst  as  far  as  local  normal  pressures  on 
the  pile  structure  are  concerned,  and  that  the  maximum  force 
can  be  as  high  as  three  or  four  times  the  frozen-in  case. 
Of  course,  the  total  lateral  load  remains  the  same. 
6.5  Locations  of  Maximum  Stress  and  Strain 
As  a  final  use  of  the  model,  information  from  the  loading  of 
Figure  6-6   will  be  presented.  Physically,  this  edge  loading 
might  be  produced  by  a  wave  driving  force  impinging  on  the 
circular  floe.  No  failure  mechanism  is  included,  since 
some  investigators  feel  that  failure  is  ductile,  caused  by  a 
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FIGURE  6-8 

Comparison  of  normal  (radial)  stress  at  pile  -  ice  sheet 
interface,  with  and  without  friction.  Loading  of  Figure 
6-6  applies. 
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limiting  strain,  and  no  reliable  theory  is  available. 
Therefore  no  estimation  of  a  maximum  total  force  will  be 
offered,  but  a  good  representation  of  stress  and  strain 
distribution  can  be  presented. 

Figure  6-9  depicts  the  variation  of  normal  (radial) 
and  shear  stresses  throughout  the  ice  sheet  and  at  selected 
angles  G  .   The  greatest  radial  stress  is  at  9=0°  as 
would  be  expected.   The  crossover  to  tension  occurs  at  6=90°. 
The  maximum  shear  stress  occurs  at  about  5^°,  again  at  the 
inner  face. 

Figures  6-10  and  6-11  show  the  variation  in  vertical  and 
hoop  strain  around  the  pile.   If  tension  strain  failure  is 
the  failure  mechanism,  it  can  be  seen  that  failure  will  occur 
in  the  sheet  vertically  (a  horizontal  crack)  and  horizontally 
(vertical  crack)  at  angles  of  0°  and  5^°  respectively.   Shear 
strain  follows  the  shear  stress  and  in  this  case  is  2/3  the 
O^e  value.   Shear  strain  is  shown  in  Figure  6-12.   The  values 
of  these  strains  are  plotted  in  Figures  6-10,  6-11,  6-12  only 
at  the  interface,  since  this  is  the  location  of  their  maxima. 
It  should  be  pointed  out  that  the  basic  shape  is  maintained 
some  way  into  the  ice  sheet,  so  that  a  crack  vrould  be  init- 
iated at  the  ice  -  pile  interface  and  propagate  radially 
outward. 

It  would  seem  heuristically  that  the  cracks  at  5^-° 
would  occur  first  due  to  the  greater  magnitude  of  the  strain. 
It  is  then  possible  that  the  crack  at  0°  will  then  follow  due 

the  redistribution  of  stress.  Both  types  have  been  reported  [15]. 
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FIGURE  6-9 

Elastic  behavior  -  Loading  of  Figure  6-6  -  Variation 
of  normal  (radial)  and  shear  stresses  throughout  the 
ice  plate 
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FIGURE  6-10 

Elastic  behavior  -  maximum  vertical  strain  location  on  ice  —  pile 
interface.   Plot  of  vertical  strain  vs.  angular  position  at  radius 
equal  to  1.  Loading  of  Figure  6-6  applies. 
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FIGURE  6-11 

Elastic  behavior  -  maximum  hoop  strain  location  on  ice  -  pile 
interface.  Plot  of  hoop  strain  vs.  angular  position  at  radius 
equal  to  1.  Loading  of  Figure  6-6  applies. 
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FIGURE  6-12 

Elastic  behavior  -  maximum  shear  strain  location  on  ice  -  pile 
interface.   Plot  of  shear  strain  vs.  angular  position  at  radius 
equal  to  1.  Loading  of  Figure  6-6  applies. 
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CHAPTER  7 
Conclusions  and  Recommendations 

7.1  Conclusions 

The  model  developed  in  the  preceeding  chapters  is 
applicable  to  a  circular  elastic  or  viscoelastic  sheet 
surrounding  a  rigid  circular  inclusion.  The  conditions  of 
plane  stress  are  assumed,  although  plane  strain  could  "be 
very  easily  implemented.   Loading  can  be  by  displacement 
or  stress  at  either  boundary,  although  mixed  conditions 
at  a  particular  boundary  are  not  allowed.   The  model  is 
favorably  disposed  for  inclusion  of  failure  criteria,  but 
none  are  included  or  proposed.   Tension  is  required  at  the 
interface  to  support  a  valid  operation  and  this  corresponds 
to  a  frozen-in  situation.   This  model  can  accept  time  dep- 
endent loadings  below  dynamic  significance. 

The  model  is  very  thrifty  with  computational  resources 
and  produces  a  very  high  degree  of  correlation  with  existing 
analytical  solutions.   The  model  is  capable  of  considerable 
expansion. 

Some  results  with  a  fictional  loading  predict  frictional 
effects  are  present  only  at  a  small  region  immediately  pre- 
ceding the  pile.   If  angular  movement  is  allowed,  the  resul- 
tant loss  of  shear  requires  an  increase  in  normal  stresses. 
If  a  separation  of  boundary  is  to  occur,  and  occurs  where 
tension  cannot  be  supported,  then  the  contact  angle  will  be 
1  90°. 
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The  range  of  time  where  ice  can  be  considered  elastic 
(below  yield  levels)  with  a  uniaxial  load  is: 

time  <  f-^1 

Other  loadings  could  be  analyzed.   Rate  of  loading  is 
important  insofar  as  time  to  failure  and  required  force. 
Finally,  viscoelastic  behavior  predicts  ice  failure  (and 
maximum  forces)  for  loading  well  below  immediate  fracture. 

Lastly,  it  is  possible  to  qualitatively  see  failure 
by  limiting  strain  or  shear  at  locations  in  the  ice  sheet 
that  are  physically  relevant.  Without  knowing  the  values 
of  biaxial  strengths  of  the  ice  material,  it  is  not  possible 
to  predict  a  maximum  force  before  fracture,  but  it  is  likely 
that  a  limiting  strain  or  shear  strength  precedes  ultimate 
crush  failure. 
7.2  Recommendations  for  Future  Work 

The  model  developed  in  this  thesis  is  extremely  rudi- 
mentary and  can  be  used  for  pile  -  ice  sheet  interaction 
under  only  the  most  generous  conditions.   Improvements  can 
be  made,  however,  to  increase  the  usefulness  into  problems 
with  more  physical  application.   These  improvements  by  order 
of  importance  are: 

(1)  Driving  Forces 
Ice  sheets  and  ice  floes  are  driven  by  the  effects  of  wind 
and  current  shears.   It  would  be  desireable  to  develope  an 
analytical  or  numerical  procedure  whereby  these  shears  can 
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"be  approximated  into  edge  loading  conditions.   The  size  of 
the  ice  sheet  considered  in  this  model  would  require  negli- 
gible wind  or  current  shear  itself,  so  that  the  two-dimensional 
approximation  could  remain  valid. 

(2)  Failure 

The  failure  mechanism  or  mechanisms  can  be  easily  applied  to 
the  model  if  a  good  choice  could  be  found.   To  date,  inves- 
tigators have  proposed  Tresca  or  Von  Mises  type  criteria, 
failure  by  limiting  strain,  rate  dependent  crush  strength, 
among  others.  The  behavior  of  ice  is  poorly  understood  in 
biaxial  loading  and  experimental  work  needs  to  be  carried 
out  to  properly  extend  the  model  into  realistic  and  workable 
significance. 

(3)  Two-dimensional  Extension 

The  model  should  be  given  another  co-ordinate  in  the  plane 
such  that  more  general  geometries  can  be  analyzed.   Of  course 
this  could  mean  larger  quantities  of  data  to  be  handled  and 
an  attendant  loss  of  simplicity.   The  problem  of  the  mixed 
boundary  conditions  and  different  shapes  of  the  ice  floes  or 
structures  could  be  adjusted  to  fit  physical  requirements, 
with  a  high  degree  of  confidence. 

(H-)     Multiple  Structure 
The  model,  once  in  a  more  general  two-dimensional  state, 
could  be  extended  to  study  the  effect  of  multiple  structures 
such  as  support  pilings  on  a  pier  or  multilegged  towers. 
Possible  reductions  in  forces  by  interference  could  be  shown. 
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(5)  Other  Material  Behavior 
Once  in  a  two-dimensional  form,  the  model  could  be  extended 
into  non-linear  and  non-elastic  behavior  in  a  straightforward 
manner.   The  one -dimensional  form  presented  in  this  thesis 
is  valid  only  where  superposition  is  valid. 

(6)  Other  Structural  Behavior 
Experimental  work  with  the  model  could  be  extended  into 
non-rigid  analysis  for  study  of  structural  responses  and 
forced  vibrations. 

Data  on  the  material  properties  of  ice  in  a  biaxial 
state  of  stress  need  to  be  acquired   for  realistic  imple- 
mentation in  the  model.   Once  this  and  the  six  stated  improve- 
ments of  the  model  are  implemented,  the  behavior  of  the  ice 
sheet-structure  problem  can  be  fully  analyzed  and  more  general 
design  equations  or  curves  offered.  Additionally,  the 
influence  of  other  important  parameters  such  as  temperature 
profile  or  variation  of  ice  thickness  can  then  be  included 
or  discarded. 
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APPENDIX  A-l 
Program  Variables 


Variables  in 

Common 

Purpose 

n 

number  of  radial  spaces 

nn 

number  of  circumferential  spaces 

nnpls 

nn+1 

ind 

indicator  of  boundary  conditions 
used  (see  listing  under  load- 
My) 

m 

wave  number 

ru 

poisson  ratio 

el 

Young's  modulus 

Pi 

3.1415... 

alpha 

viscous  coefficient 

rades 

number  of  wave  numbers  to  be 
processed 

rout 

outer  radius 

rin 

inner  radius 

idsp 

4xN-r  number  of  radial  display 
intervals 

ibar 

termination  variable 

Variables  Consistently  Used 

irowmax 


number  of  simultaneous  equations 
to  be  solved 
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tstop  stop  time  limit  (visco  elastic 

solution) 
tstep  time  step  interval 

time  time  since  start  (elastic:  time  =  0) 

disint  time  display  interval  of  viscoelastic 

solutions 
streslim  elastic  limit  of  stress 

strnlim  limiting  strain 

rad  radius 

hoop  hoop  stress 
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APPENDIX  A-2 
Description  of  Arrays 

Array  Subroutine 

/blokl/a (1601,17)  elastic,  zerabx,  loada,  strain, 

double  precision  loadbdy,  viscp,  dpband 

purpose:  holds  coefficient  matrix  in  band  compressed 
form  -  after  call  to  dpband  holds  the  factored  matrix 

/blok2/b_(l601,l)  elastic,  zerabx,  loadbdy, 

visco,  loade,  boundb,  strain, 
dpband 
purpose:  solution  vector  of  the  system  of  linear 
equations  represented  by  a 

/blok3/c (204,21)  elastic,  zercdr,  recomb, 

writr2 
purpose:  stores  total  recombined  solutions  at  every 
eighth  node  (row)  and  all  circumferential  positions 

/blok4/d(ll,8)  elastic,  zercdr,  fourier, 

loadbdy,  writrO 
purpose:  stores  fourier  coefficients  of  boundary 
conditions  for  0  through  10  wave  numbers  (row),  inner 
conditions  (first  four  columns),  outer  conditions 
(second  four  columns)  as  determined  by  variable  ind 

rl(8,21)  loader,  fourier 

purpose:  stores  boundary  information  from  main  before 
fourier  decomposition 
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ws(21)  fourier,  simps 

purpose:  working  vector  to  be  integrated  by  simps 

/blok6/xl (14436)  elastic,  zerabx,  dpband ,  visco 

double  precision 

purpose:  working  storage  used  by  dpband 

/blok8/r(204,ll)  elastic,  zercdr,  recomb 



purpose:  stores  elastic  wave  number  solutions  for  every 

eighth  node 
/blok9/t(51*H)  elastic, zercdr ,  recomb 

purpose:  stores  hoop  stress  wave  number  solutions  for 

every  eighth  node  (similar  to  /blok8/) 
/bloklO/tt ( 51 , 21 )  elastic,  zercdr,  recomb,  writr2 

purpose:  stores  total  recombined  elastic  hoop  stress 

solution  (similar  to  /blok3/) 
/blokll/rr(l6o4,l)  visco,  loade 

double  precision 

purpose:  working  vector 
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APPENDIX  A -4 
User  Notes 

Program  Operation 

The  program  flow  is  described  in  Figures  4-10  and  5-2. 
The  material  constants  E, H ,  c*  are  set  in  either  the  main 
calling  program  or  subroutine  elastic  prior  to  compilation. 
In  the  listing  presented  in  A-3,  E  and  o<  are  non-dimensional 
with  a  value  of  1.0.   The  program  is  not  restricted  to 
non-dimensional  solution,  however,  and  any  consistent  values 
could  be  used. 

User  inputs  are  via  keyboard  and  consist  of  3  inputs. 
For  elastic: 

1)  mdes  -  after  reviewing  fourier  coefficients  of 
boundary  conditions,  the  user  is  asked  for  the 
number  of  fourier  coefficients  to  be  used  in  the 
solution.   This  feature  is  incorporated  to  save 
calculation  of  meaningless  higher  wave  number 
solutions. 

2)  The  user  is  queried  on  whether  a  viscoelastic 
solution  is  required.   This  feature  is  used  to 
permit  an  elastic  solution  only. 

3)  Display  of  theta  information  -  This  feature 
automatically  combines  elastic  wave  solutions  and 
calculates  circumferential  variation  for  intervals 
of  .05  TT  and  radial  variation  of  50  intervals. 


142 


Display  here  is  set  for  .01 TT  and  radial  variation 
o£  ten  intervals. 
If  a  viscoelastic  solution  is  desired,  additional 
keyboard  inputs  are : 

1)  tstop  -  desired  stop  time,  entered  in  units  of 
non-dimensional  time. 

2)  tstep  -  interval  calculation  time  -  This  is 
important  as  a  small  step  must  be  chosen  to  ensure 
convergence. 

3)  disint  -  time  display  interval  of  wave  number 
solutions  -  It  is  anticipated  that  not  all  small 
discrete  time  intervals  will  require  display. 
disint  should  be  an  integral  multiple  of  tstep. 

Loading  Input  -  rout  must  be  set  in  the  main  program.  The 

vectors  putl  and  put4  contain  boundary  condition  loading 
information  as  determined  by  ind,  also  ^ent  from  main. 
putl  -  put4  contain  diecrete  data  for  loading 
circumferentially  through   0=TT  .   (The  other  half 
is  identical  by  symmetry). 

Inner  Boundary  Conditions 


ind 

putl 

put2 

1 

Uy. 

u© 

2 

Or 

6r© 

3 

U-r 

u© 

4 

<% 

6r© 

Outer 

Boundary  Conditions 

put3 

put4 

Uy 

^© 

6v 

e>v» 

6, 

Gte 

u* 

ue 

Caution  must  be  exercised  with  ind  =  2,  since  a  rigid  body 
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translation  may  be  produced. 

Output  of  Program  -  Three  distinct  listings  are  produced: 

1)  Fourier  coefficients  of  boundary  conditions  -  all 
wave  numbers 

2)  Wave  number  solutions  of  C>x  >§re  ,^v  »  ^s,  §e  ,  §cff 
6eff  *  node,  radius  for  upto  4-00  nodes.   In  the 
listing  of  Appendix  A -3,  only  ten  radial  intervals 
are  displayed,   idsp  will  vary  display  points. 
Wave  number  solutions  are  displayed  for  the  elastic 
case  and  at  each  time  display  intervals  determined 
by  disint.   The  headings  "stress"  and  "strain" 

contain  an  effective  stress  and  strain  as  defined 
by:      ,  i 


€-c^*fe  *  ms  -**;)1 


3)  Stresses  by  angular  position  -  contains  up  to 
fifty  radial  intervals  with  the  total  recombined 
fourier  series  representation  of  the  elastic 
solution  displayed  at  intervals  of  .OlTT  . 
Other  Notes 

Viscoelastic  boundary  conditions  (rates)  are  preset 
to  zero  in  this  program,  but  are  not  necessarily  restricted 
to  this  value  or  even  a  constant  value.   Rate  information 
must  be  in  size  with  integration  of  time  step  and  convergence 
must  be  ensured  overall. 
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streslim  and  strainlim  are  preset  constants  corresponding 
to  an  elastic  limit  and  a  strain  failure  limit.   The  solution 
is  deemed  inappropriate  beyond  these  levels  and  a  solution 
is  terminated  for  that  particular  wave  number. 

The  accuracy  of  the  fourier  subroutine  is  dependent  on 
the  circumferential  interval.  For  complex  loading  it  is 
conceivable  that  this  interval  can  be  decreased  by  increasing 
nn  and  of  course  pertinent  array  sizes. 

Storage  of  solutions  in  the  B,  C,  and  R  arrays  for  a 
particular  variable  are  in  every  fourth  row,  e.g.; 

solution 

2,mb(u),b(i+l,l),b(u2,i)1b(u3,l) 
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APPENDIX  B 

Analytic  Elastic  Solution  for  the  Case  of  a  Circular  Ring 

Wavenumber  m=0 

The  stress  function  is  known  to  satisfy  the  biharmonic  relation 

I  VA<t-*0  (B.l) 

In  a  two-dimensional  elastic  material.  For  the  case  of  no 

circumferential  variation  (as  in  wavenumber  0),  (B.l)  reduces 

4H   +l£S..^-lH ■       l.Si  (B.2) 


Let 

V  \°tr 

Then 

41-  il  £il  -  ill 

dr  "  d$   dv          V   d^[ 

dV      d  /d*r\    1    /dH    dt 

d-r^Pla-Y3    3ata4ZdiJ 

d<+     imH      fcdH.,,dH     /0a_f\ 
dT*  "-^IdJ^  "  *3T»     l  dp        Jf) 

(B.2)  becomes 

The  general  solution  of  which  is 
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Finally      ^-_  C^Mo^v  4-  Cz  lo^v  +  Cz  v1  *  C3  I03*  v  C4 

l  d4 


Now 


V© 
Strain  relations 


6v  =  t  tb  =  c*  0  +  *Aoy)  +ic2  ♦  ^2      (b.3) 
0,e-  o 


do  .  U 


elv         V 


vy  ~    elv-  ce        r  tre       3£       v 

In  plane  stress 

62--  o 

I  *-fU-*0  (B'5) 

t- •  ±  ( &•  -  *  6*)  (b.6.) 

Integrating   (B.5) 

Eu~-    ^(B.3)^    -^i(B.4)ar 

Into  £„.-.  C^O""^)  *■  1(^^)1^  \e>^v-v-)]  +  2Ci(\->J)v. 

-C3(\+^T    +  C5  (B.7) 

I  (B.6)  becomes 

I  Ea--  C^v  (!-*»)  +  2v(l-*)loyr]  V  2Cz(\-pV 

I  -C,Cu^T  (B.8) 

;For   (B.7)  to  have  the  same  representation  as    (B.8) 

I  cx-c*--  o 

Therefore,   for  wavenumber  0 

Eu  -    2  Ci  0-*0r-   C^Cl+v^—  (B.9) 


Wavenumber  m=l  154 

From  the  general  solution  for  the  stress  function  in  circular 
coordinates  presented  by  Timoshenko  [29],  that  portion  cor- 
responding to  the  first  wavenumber  is: 

*P  a   2    r  ©  Su*  ©  +  (,  b,  V3  +  a,  v"  ■*.  V>,  v-loy  )  Co*  Q 


Since 


A         1    If  '      ^1 


3v~ 

Then  (5V  *  L.U,+b,  )*"'  +  2.b,r  -  2ck!r"3]  cos© 

6©  "  ( 6\&,v-  4-2a.,v-   t^v"  )co$© 

6r6  *  (  2 b,  ir  -  2a,  v"3  *  W*  v"'  )  sC,  9 

Employing  stress-strain  (4.11)  and  integrating  with  respect  to  r 

+[io'lo^y]  ^  Cos©  V  £(©)  (B.ll) 

Employing  (4.12)  and  integrating  with  respect  to  O 

-  f^C©)el©  +  FO)  (B.12) 

Prom  (Boll)  T   3©"  =  '  Ua'+  b'  J  r    +  V  *  *i «"  JSu*e  +   , 

(B.13) 

Prom  (B.12)   C  f^  •  [-  (^^  +  >°^  -  2a,  ^]  *fc  ft  + 

(B.14) 
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And  also  from  (B.12) 

+  r  jcgijQ  _  _l  (B.15) 

)      v  y. 

Equating  (4,8)  and  (4.5)  and  using  (B.13),  (B.14),  and  (B.15) 
or         ^(e)*  ^ce)ae  *(v*'-f).  [41b,,*(»-^0h]iU»O 

(B.16) 
(vF-F1)  must  necessarily  be  zero 

V-F*-  F  =  O  =>  F=  Ct- 
Differentiating  with  respect  to   © 

fM  •  [^  (l-^a,]cos©         (B#17) 
The  solution  to  the  differential  equation  (B.17) 


Wavenumber  >  2  ,- 

From   [29] 


+ 


+  (w-iVva-i)  bv>  v*"*^  cos*  9 

-v\(v\-0lo*  /  ]s6*te© 

Again  employing    (4.11)   and   integrating 

t  Uy  ■    -  Lv*A«»,r     +^-2Jl°M*>      -IaA^V        -(nvljkv.v       J* 

(B.18) 
Employing    (4.12)   and  integrating 

.  siaa  w€>  J-  v->  Lu  A*  v      ^^r       ^fl.i»V 

(B.19) 
Prom   (B.18)       £  ^ii 


-w-l 


-V-     I  +y*  )>  a.*  r      J-  *(>vlj  v>v,v  -  v\  <X*  v 


(B.20) 


Prom  (B.19)       E  ^  »  f[wM^v^+  C^'X^bvA  vi(u*iW  v"^ 

.  r*  _  w ( U4 f) ftj  v"U"l  \a (>-0 *o h *** ][c S0u w ©  vp 

(B.21) 
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Also  from   (B.19) 

+  *v>,w*"*]wwe  +  ^  ^tcftd*-£.  (B.22) 

Equating   (4.8)   and    (4.9)   and  using   (B.20),    (B.21),   &   (B.22) 

v 

r  •  cvOt 

■£   r      A  Si~  ©    +     B  6**   0 


These  general  solutions  must  "be  employed  with  the  correct 
"boundary  conditions  to  find  the  unique  constants  associated 
with  each  wavenumber. 
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APPENDIX  C 
Buckling  of  a  Circular  Ring  Subject  to  Stresses  Uniformly 

Distributed  Around  the  Edge 
The  loading  applied  to  the  condition  depicted  in  Figure 
C-l  is  analyzed  for  the  first  buckling  mode.  Assuming  the 
deflection  surface  is  a  surface  of  revolution,  Timoshenko 
[28]  gives  the  required  differential  equation: 

g±  d*  _    _   Qv^  (C.l) 

r  dyl     it-    T     D 

where  Q  is  the  shearing  force/length 

and  D  is  the  flexural  rigidity 

Q«N*sW4>  =  N„4>  (c-2) 

Let         !!!*  *  o<Z  (C.3) 

D 

(C.3)  becomes 

**££♦»£  ♦<-V-0*«°       (o.») 

If  u--  <*v-  (variable  change) 

Then        ^^l  +  u^  +  C^-O-O         (o.5) 

The  general  solution  of  which  is 

<b«  A.JJCuVAaY.C*)  (c.6) 

Applying  boundary  conditions 

^>(v-»<x)  =  o     clamped  inner  edge      (C.7) 
f di  4.  p  ^.V  O   v-*k  free  outer  edge   (C.8) 
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Figure  C-l 


Buckling  of  a  thin  circular  plate  around  a  rigid 
inclusion 
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sinCe         ^.*  A,c*I,'(u)  4-  Az*  Y '(".")  (c-9) 

The  boundary  conditions  (C.7)*  (C.8)  become 

A,J,Ucx)  +  A1Y1Co<aO  =  o  (cio) 

*(A,J/U^A2Y/(^))  + 
+  ^(A,Tl(-b)+Ai.Yluy).o        (c-11) 

(cio)  is         Ai«-A,    '        x  (c-12) 

And   (C.ll)  <X(J,U)-  ~^V'(^))  +  -£(T,(^)-    (C13) 

Y,  (**)    '  y 

(C.12),    (C.13),    (0.14)    are  combined   into  a   single   equation. 
The  first  zero  appears  at  about  a  value  of   .42 

To  apply  aa,>      h   z  20    ■>  ^'"3 

Where  U-.4Z/0-     then         Nt*~.n&D 

Q  =     ^  ,   the   f lexural  rigidity 


Yields  6cr  =  —    S  -0lfoS  EW1/*1 

In 
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